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Summary 


The influence of thermal nonequilibrium in rarefied hypersonic 
flows has been considered. In a compressed flow, such as a shock 
wave, typical ^relaxation effects cause the rotational and vibra- 
tional temperatures of molecules to lag behind that for the trans- 
lational energy mode. In the energetic flowfields which surround 
reentering spacecraft thermal nonequilibrium effects are most im- 
portant due to their influence on the chemical activity of the gas. 
The application of Computational Fluid Dynamics techniques to such 
flowfields must be treated with great care. At altitudes greater 
than about 80km above the Earth' s surface the rarefied nature of 
the ambient atmosphere makes the use of traditional continuum tech- 
niques inappropriate and resort must be made to particle simula- 
tion methods , The manner in which energy exchange between the var- 
ious energy modes is modelled in the particle simulation is clearly 
of great importance . A new model for calculating the probability 
of energy exchange between the translational and rotational modes 
has been developed. A full description of this work is contained in 
Appendix A which has been submitted for publication in The Physics 
of Fluids. The role of vibrational energy exchange has also been 
investigated. A model for determining the probability of energy 
transfer to the vibrational modes has also been derived. This model, 
together with that described in Appendix A, have been incorporated 
into a one-dimensional study of the flow of air past a blunt body 
at an altitude of 90km. The results of this study are discussed in 
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Appendix B. This material was presented at the 20th AIAA Fluid Dy- 
namics, Plasma Dynamics, and Lasers Conference, Buffalo, NY, June 
1989, and has been submitted for publication to The Journal of Ther- 
mophysics and Heat Transfer. 
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APPENDIX A 



ROTATIONAL-TRANSLATIONAL ENERGY TRANSFER 
IN RAREFIED NONEQUILIBRIUM FLOWS 


Iain D. Boyd 

Eloret Institute, NASA Ames Research Center 
Moffett Field, CA 94035 

Abstract 

A new model for simulating th.e transfer of energy between the trans- 
lational and rotational modes is derived for a homogeneous gas of di- 
atomic molecules. The model ha£ developed specifically for use in 
discrete particle simulation methods where molecular motion and inter- 
molecular collisions are treated at the molecular level. In such methods it 
is normal to assume a constant rotational collision number for the entire 
flowfield. The new model differs in that a temperature dependence is in- 
troduced which has been predicted by theory and observed in experiment. 
The new model is applied to the relaxation of rotational temperature, and 
is found to produce significant differences in comparison with the model 
normally employed at both high and low temperatures. Calculations have 
also been performed for a Mach 7 normal shock wave. Large differences 
in the solutions are again observed, with the new model offering an im- 
proved correspondence to the available experimental data. 


PACS Numbers 

34.50.Ez Rotational Energy Transfer. 
47.45. -n Rarefied Gas Dynamics. 
02.70.+d Computational techniques. 
47.40.Nm Shock waves. 
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I. Introduction 


The use of discrete particle simulation methods in fluid mechanics 
has become widespread with the recent increase in interest in rarefied, 
nonequilibrium gas flows. This interest has mainly been concerned with 
high altitude rocket and thruster expansion plumes 1,2 , and with the flow- 
field surrounding a space vehicle as it re-enters the Earth’s atmosphere at 
high velocity 3 . The most commonly applied particle simulation tech- 
nique is that developed by Bird 4 and is termed the Direct Simulation 
Monte Carlo method (DSMC). This solution scheme proceeds by mod- 
elling all physical and chemical phenomena at the molecular level, and 
is therefore preferred to more traditional continuum formulations when 
rarefied effects become important. 

One of the most significant effects observed in rarefied flowfields 
is the presence of a large degree of nonequilibrium between the various 
internal energy modes of the gas. It is therefore of great importance that 
such phenomena are adequately simulated within the DSMC technique. 
Presently, the exchange of energy between translational, rotational, and 
vibrational modes in the DSMC method is treated in a simplistic manner 
which relies more on convenience than any physical basis. In particular, 
the existing models fail to allow the collision numbers associated with the 
rotational and vibrational modes to exhibit any temperature dependence. 

The purpose of the present work is to improve the modelling of 
translational-rotational energy exchange for use in discrete particle simu- 
lation methods. A temperature dependent collision number employed in 
continuum calculations is used as a basis to develop a model which re- 
produces equilibrium behaviour in results obtained with the DSMC tech- 
nique. The new model is then applied to a simple relaxation problem, and 
to a strong shock wave in which a significant temperature range exists. 

II. Translational-Rotational Energy Exchange 

In a gas of colliding molecules, energy is continually being trans- 
ferred between the various internal modes. These collisions tend to push 
the internal energy distributions towards their equilibrium state. The 
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number of collisions required to equilibriate a particular mode of each 
molecule is called the collision number Z of that mode. Each internal 
mode has a separate collision number, and it is generally the case that 


E translation < ' Z rotation < ' ^vibration 


It therefore takes a greater number of collisions for the rotational 
mode to reach equilibrium than the number required for translational en- 
ergy. Thus, if the collision rate in any reference volume in the flowfield 
does not allow each molecule to collide Z rotat i on times before leaving, 
then nonequilibrium exists in that volume. In the treatment of rotational 
energy exchange between the translational and rotational modes, the fol- 
lowing expression is often employed in continuum calculations: 


dEft E*p — E R 

~r = — — - (i) 

dt Tr 

where E R is the rotational energy, E* R is the equilibrium value, and r R 
is the rotational relaxation time which is generally regarded as being a 
function of temperature and pressure (or density). 

While Eq.(l) is suitable for continuum calculations, it is inappropri- 
ate for use in particle simulation schemes in which all macroscopic flow 
properties are derived by averaging molecular properties over large sam- 
ple sizes. The particle simulation technique considered here is the Direct 
Simulation Monte Carlo method. In this technique, the large number of 
molecules in a real gas are simulated by a much smaller representative 
set of model particles. The motion of these particles is followed through 
physical space and intermolecular collisions are treated on a probabilistic 
basis. Several different collision schemes have been developed and most 
of these have been categorised by Nanbu 5 . The scheme employed in the 
present work is the Time Counter method of Bird which has the com- 
bined advantages of numerical efficiency and low statistical fluctuations 
(see Ref. 6). 

In addition to retaining molecular velocities, it is possible to sim- 
ulate the internal energy modes associated with polyatomic gases. In 
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order to avoid the complexities of computing individual transitions be- 
tween the various internal quantum states, a continuous distribution of 
internal energies is usually assumed. While the rotational energy of a 
molecule should be simulated by considering its angular velocity vector, 
it is computationally simpler to ignore all effects due to the orientation 
of the colliding bodies. Thus, it is usual to associate a scalar rotational 
energy value with each molecule. 

Rotational energy exchange schemes involving classical laws of me- 
chanics have been proposed by MacPherson 7 and individual quantum 
level transition probabilities for nitrogen have been calculated by Deiwert 
and Yoshikawa 8 . Both of these methods suffer from computational over- 
heads which prohibit their use in discrete particle simulations of complex 
flowfields. 

In a particle simulation scheme it is more economical to employ the 
concept of a collision number in preference to either Eq.(l) or the work 
described in Refs. 7 and 8. The rotational collision number is usually 
defined as 

Zr=— ( 2 ) 

T c 

where r c is the mean time between collisions. The transfer of energy be- 
tween the various modes is usually implemented into the DSMC tech- 
nique using the Larsen-Borgnakke phenomenological model 9 . In this 
model, a fraction 0 of all collisions are considered to be inelastic and 
the remaining (1-0) collisions are treated as perfectly elastic spheres. If 
the collision number Z of an internal mode is given by Eq.(2), then the 
probability of energy exchange involving this mode, 0, is conveniently 
defined as 


+ m f (3) 

where Z t is the translational collision number. It is usual in DSMC calcu- 
lations to take Z t to be equal to unity and Z R = 5, so that a constant value of 
(f) R = 0.2 is obtained. This value for Z R has been chosen to coincide with 
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results obtained from ultrasonic 10 and shock wave 11 experiments. The 
calculation of an inelastic collision is then performed by statistically sam- 
pling values using the equilibrium distribution functions together with the 
total collision energy (including internal mode contributions). 

- By assuming that the rotational collision number is constant over the 
entire flowfield, any temperature dependence is thereby neglected. This 
concept is in direct opposition to the theoretical findings of Parker 12 , and 
Lordi and Mates 13 , and the experimental data of Camevale et al 14 . Parker 
has used an empirical non-impulsive potential model which incorporates 
a small degree of asymmetry to derive an expression for the rotational 
relaxation time t r . By initially assuming that there is zero energy in 
the rotational modes of a gas of homonuclear diatomic molecules, the 
following approximate expression is obtained: 



where T* is the characteristic temperature of the intermolecular potential, 
and (Zfl)^ is the limiting value. 

While Parker’s expression is derived from an analysis involving a 
large number of assumptions, the temperature dependent nature of Eq.(4) 
is in agreement with the more rigorous treatment of Lordi and Mates 
who performed classical trajectory calculations. In the current work the 
value of (7'r)oo is chosen so as to obtain the best correspondence between 
Parker’s results and those of Lordi and Mates. 

A previous attempt to introduce a temperature dependent nature of 
the rotational collision number into DSMC calculations was reported by 
Davis et al 15 . It was proposed that Z R should be specified as a func- 
tion of the total collision energy of the colliding molecules, i.e. the sum 
of translational and rotational collision energies. The exact derivation 
of the expression employed is omitted and the model fails to satisfy the 
principle of detailed balance. It is therefore the intention here to derive 
an expression from Eq.(4) which may be utilized in discrete particle sim- 
ulations. 
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III. Probability of Translational-Rotational Energy Exchange 

In the following, the principal aim is to derive an expression for the 
reciprocal of the rotational collision number from the molecular collision 
properties of the gas. A decision must therefore be made on the collision 
parameter best suited for this purpose. As Eq.(4) is strictly a function of 
translational temperature this suggests that the relative velocity of colli- 
sion, c r , should be employed. For an equilibrium gas, it may be shown 
that the average value of a quantity Q which is a function of c r alone may 
be expressed as 


< Q >= 


? f-UhA 2-i 
_ \2kTJ_ 

r(2 — w) 



Qc r 3 2u} exp 


—m r c T 2 \ 
2kT ) dCr 


(5) 


where m r is the reduced mass of the collision and w is a constant used in 
the Variable Hard Sphere (VHS) collision model of Bird 16 and is related 
to the viscosity temperature exponent Q by u = £ — 0 .5 . In the VHS 
model the collision mechanics are treated in the same manner as for hard 
spheres while the total collision cross section is a function of the relative 
velocity and is given as 


(Jj. = 


( 2 j 

1^(2 - W )^r 0 J 


( 6 ) 


where cr 0 is a reference collision cross section defined at temperature T 0 . 
In this formulation, a value of u> = 0 corresponds to the hard sphere inter- 
molecular potential and u = 0.5 gives the Maxwell molecule. 

The probability that a molecule transfers energy between its trans- 
lational and rotational modes due to a collision is given by the reciprocal 
of Eq.(4) which may be written as 


(pR = ai + 


Q2 03 

T x i + T 


(7) 
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where ai , a 2 , a 3 are constants. To reproduce the global expression of 
Eq.(7) in a discrete particle simulation, an appropriate form of Q must be 
found for the integration shown in Eq.(5). Let us assume that each of the 
three terms on the right hand side of Eq.(7) may be obtained in this way 
be choosing a general condition Q=bc r n . Thus 


< Q > 


_ 2 b( 


T(2 - 


2 —U) /*C 

Tj_ / 

- w) Jo 


Tn r N 2 —u) poo 

2 kTl I _ n+3— 2a; 

c r 


exp 


— m r c T 
2kT 


d,Cr 


( 8 ) 


By making the substitution x = Eq.(8) becomes 


2kT) oo 

< Q > = r ™ r - 7 [ x% +l ~ u exp(—x)dx 

T( 2 -uj) J 0 

1(2 — w) \ m r / 

On examination of Eq.(9)’it becomes clear that the values of n which 
must be employed for each of the terms in Eq.(7) are 0, -1 and -2 respec- 
tively. Therefore, the full expression for the average value of probability 
for energy exchange between the translational and rotational modes is 


<f>R 


(Zr) oo 

Zt ' 


T ( 2 — lo) / 2 kT* \ 2 7T ^ 
T( | — uj) \m r c r 2 ) 2 

T( 2 - uj) ( 2 kT^_ \ (ifj_ 

F ( 1 - w) \m r c r 2 / V 4 



( 10 ) 


The form of Eq.(10) shows that there is a singularity at c r = 0. How- 
ever, the collision probability for a pair of molecules which have zero 
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relative velocity is itself zero, so that this case does not present a prob- 
lem. The application of this equation in a DSMC calculation to a number 
of test cases is considered in the following Sections. 

IV. Calculations 

To assess the validity, and the significance, of the energy exchange 
model proposed in the previous Section, the expression in Eq.(10) has 
been implemented into a number of DSMC calculations. In the present 
study, a value of Z x =1.5 has been employed in accordance with the result 
of Fritzsche and Cukrowski 17 . The first problem undertaken involves a 
heat bath of nitrogen molecules which are initially in thermal equilibrium. 
The values employed in Eq.(4) for nitrogen which best correspond to the 
calculations of Lordi and Mates are T*=91.5K and (Z R ) oo=23.5. To be- 
gin with, Eq.(10) was employed in such a way that the relative velocity 
of a given collision determined the instantaneous probability of energy 
transfer. However, such a procedure fails to preserve the equipartition of 
the translational and rotational modes. This is because energy exchange 
then occurs for collisions which have a lower than average translational 
energy. It was therefore necessary to employ a value for (j) R which is aver- 
aged over every collision regardless of whether energy exchange actually 
occurs. This procedure has been tested over a range of temperatures and 
some example calculations for the values of Z R derived from the DSMC 
results together with Eqs.(3) and (10) are shown in Fig. 1. It is clearly 
shown that the results obtained at the molecular level are found to match 
Eq.(4) almost exactly. 

Having established the fact that Eq.(10) reproduces the continuum 
result under equilibrium conditions, it is relevant to consider the appli- 
cation of the model to some nonequilibrium problems. The first of these 
models translational and rotational temperature relaxation in a homoge- 
neous gas. Initially, the rotational energy is set to zero, and the transla- 
tional temperature is given by 5/3 of the final equilibrium value. As time 
marches forward in the simulation, the temperatures in the translational 
and rotational modes will converge to the equilibrium result. This re- 
laxation problem is therefore unsteady, so that the reduction of statistical 
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fluctuations is achieved by averaging results over many complete runs 
from start to finish. 

Relaxation profiles obtained with (f> R = 0.2 and Eq.(10) are shown in 
Fig. 2 for rotational temperature. The working gas is nitrogen, and time 
is normalized with respect to the mean time between collisions at the 
equilibrium temperature which is 300K in this case. It is seen that small 
distinguishable differences occur in the results calculated by the two en- 
ergy transfer models. For the translational temperature range considered 
in this calculation, the value for <f> R obtained from Eq.(10) is greater than 
0.2. The rate of transfer of energy to the rotational modes is therefore 
greater, and the rotational temperature converges to the equilibrium value 
in a shorter time. Another important feature of the results shown in Fig. 
2 is that equipartition of the thermal modes is achieved at equilibrium 
with the new model. In Fig. 3 the rotational relaxation history is shown, 
this time for a final temperature of 5000K. Under these conditions, the 
values of <f> R calculated with Eq.(10) are much less than 0.2, and so the 
relaxation process is slower. 

Unfortunately, the results shown in Figs. 2 and 3 cannot be veri- 
fied experimentally. However,, these theoretical calculations indicate that 
significant differences occur due to the introduction of a temperature de- 
pendent collision number. It is therefore appropriate to make calcula- 
tions of a flow for which rotational temperature measurements are avail- 
able. Assessment has been made of the differences found for the solu- 
tions obtained with the two models in the calculation of a strong shock 
wave. The two different sets of conditions investigated are listed in Table 
I and the working gas is again nitrogen. The first set of conditions has 
been chosen to match the experimental data of Robben and Talbot 11 who 
produced a number of shock waves in a wind tunnel by expanding a jet 
through a shock-holder. The high Mach number shock waves studied by 
Robben and Talbot are much thinner than those produced by Alsmeyer 18 
in a shock-tube. The explanation for this phenomenon lies in differences 
in the upstream flow conditions in the two experiments. In Alsmeyer’s 
study, the upstream conditions should ensure the existence of thermal 
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equilibrium. To obtain high Mach number shock waves, it was necessary 
for Robben and Talbot to allow the expanding jet to undergo significant 
rarefaction. In such expansions the translational temperature will nor- 
mally relax more efficiently than the rotational mode (see Ref. 2). Thus, 
at the upstream boundary, the translational temperature will probably be 
less than the value of 28K obtained for the rotational temperature. 

In the present calculations, thermal equilibrium is assumed at the 
upstream boundary. To reproduce the thin shock wave of Robben and 
Talbot it is necessary to employ an unusually low viscosity temperature 
exponent. Such an undertaking leads to a reduction in the molecular colli- 
sion rate which is quantitatively consistent with the upstream translational 
temperature being lower than that for rotation. 

The simulation procedure employed follows that of Bird 19 . The pro- 
files for density and temperature are shown in Fig. 4, and have been nor- 
malized as follows: 


^ P~ P i 

P = 

P2 - pi 


( 11 ) 


T = 


T-T i 
T 2 -Ti 


( 12 ) 


The calculated density profiles obtained with each of the energy 
transfer models were found to be exactly consistent. Therefore, only 
the result obtained with Eq.(10) is shown, and it is evident that excel- 
lent agreement is found with the experimental data. The reciprocal shock 
wave thickness for the density profile is about 0.41 which is in good 
agreement with the value derived by Alsmeyer 18 for the data of Robben 
and Talbot. The correspondence between the theoretical and experimen- 
tal results for rotational temperature is not as good. Significantly, the pro- 
file obtained with the new model shows an improvement of up to 20% on 
that calculated with <f> R = 0.2. The results shown in Fig. 4 indicate that 
the new energy transfer may be successfully applied in a discrete particle 
simulation, and gives better correspondence to experimental data. 
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A further set of calculations has been made for a Mach 7 shock 
wave where the upstream temperature is taken to be 300K and less rar- 
efied conditions are assumed. This allows assessment of the effect of 
the new energy transfer model when the collision number is greater than 
5 and gives conditions similar to those investigated by Alsmeyer 18 us- 
ing a shock tube. The assumption of thermal equilibrium in the upstream 
portion of the flowfield is now more certain. The normalized density pro- 
files obtained by implementing Eq.(10) for the two sets of conditions are 
shown in Fig. 5. The effect of simulating a shock tube study is to increase 
the shock layer thickness and occurs because of differences in the rate of 
rotational energy exchange and viscosity effects. The calculated value 
for the reciprocal shock thickness again shows good agreement with that 
which may be deduced from Ref. 18. 

The normalized profiles for translational and rotational temperature 
obtained from these simulations are shown in Fig. 6. For the range of 
temperature found in this simulation it is observed that more significant 
differences exist for the two solutions. As expected, the rate of transfer 
of energy to the rotational modes is reduced by application of Eq.(10). It 
should be noted that the differences between the results obtained with the 
two energy transfer models will increase as the translational energy in the 
flowfield is increased. It is therefore proposed that the new model may 
have significant impact on re-entry calculations such as those presented 
in Ref. 3 where the translational temperature may exceed 30,000K. 

Concluding Remarks 

The introduction of a temperature dependent collision number into 
discrete particle simulation schemes has been found to be significant. 
Even for the simple case of relaxation in a homogeneous gas differences 
are discernible between the results obtained with the new model and those 
calculated under the assumption that the collision number is everywhere 
equal to 5. The application of the new model to a shock wave problem un- 
der very rarefied conditions gives improved correspondence with experi- 
mental data in comparison to a constant collision number. It is concluded 
that the new model should be preferred in calculations which involve 
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a significant degree of rotational nonequilibrium. Finally, it should be 
noted that the model presented here may also be used to allow separate 
rates of translational-rotational energy exchange between the different 
molecules in a gas mixture. 
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Table I. Upstream conditions for Mach 7 shock wave calculations. 


pi (kg m 3 ) 

Ti (K) 

ui (ms -1 ) 

3.12xl0- 5 

29 

764 

7.77xl0 -3 

300 

2470 


Figure Legends 

Fig. 1 Rotational collision number obtained from Eq.(3) and Eq.(10). 

Fig. 2 Homogeneous relaxation of the rotational temperature of nitrogen 
at 300K. 

Fig. 3 Homogeneous relaxation of the rotational temperature of nitrogen 
at 5000K. 

Fig. 4 Normalized profiles of a Mach 7 shock wave. Conditions chosen 
to match those of Robben & Talbot 11 , see Table I. 

Fig. 5 Normalized density profiles of a Mach 7 shock wave for the two 
sets of upstream conditions investigated, see Table I. 

Fig. 6 Normalized temperature profiles of a Mach 7 shock wave. Up- 
stream conditions: T\ =300K, p\ =7.77xl0 -3 kg m -3 , see Table I. 
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